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We consider the problem of screening of an electrically charged impurity in a clean 
graphene sheet. When electron-electron interactions are neglected, the screening 
charge has a sign opposite to that of the impurity, and is localized near the impurity. 
Interactions between electrons smear out the induced charge density to give a large- 
distance tail that follows approximately, but not exactly, an behavior and with 
a sign which is the same as that of the impurity. 
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I. INTRODUCTION. 

With the recent explosion of interest in graphene, there are numerous experimental mo- 
tivations for understanding the influence of impurities on its electronic and transport prop- 
erties. For non-interacting electrons, the influence of a dilute concentration of impurities on 
transport properties has been investigated in some depth [1]. Here we shall instead study 
in some detail the physics associated with a single impurity carrying electrical charge Z. 
Nanoscale studies of the electronic properties of a single graphene sheet have recently be- 
come possible [2, 3], and so it should eventually be possible to observe the variation in the 
charge density and the local density of states as a function of distance from the impurity. We 
shall show here that this spatial structure is a sensitive probe of the strong correlations be- 
tween the electrons in graphene, and of the unusual nature of screening in a two-dimensional 
semi-metal with a Dirac dispersion spectrum. 

For non-interacting electrons, the influence of a Coulomb impurity exerting a potential 
Ze^/(47reor) (where r is the distance from the impurity) was studied some time ago [4]. 
This case is equivalent to the familiar "Friedel problem" but for Dirac fermions. However, 
even for this seemingly simple case, there are subtleties which were overlooked in the initial 
treatment [4], and corrected in Ref. 5. A number of papers appeared [6-8] while our paper 
was being written, presenting additional results on this non-interacting problem. We shall 
review and extend the results of Ref. 5 for non-interacting electrons in Section H. We shall 
then proceed to the full treatment of the impurity problem, and allow for electron-electron 
Coulomb interactions. 

In short, our results are as follows. For noninteracting electrons, the screening charge is 
a local delta-function in space to all orders in perturbation theory over the impurity charge. 
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The sign of this screening charge is opposite to that of the impurity, as is usually the case. 
However, once interaction between electrons is turned on, the screening charge develops 
a long-range tail, even for small impurity charges. The tail follows approximately an 
law, with a coefficient which varies quite slowly with r. Notably, the sign of this tail is 
the same as that of the impurity. The long-range tail of the screening charge, thus, is a 
sensitive probe of the interaction between electrons, in particular to the renormalization of 
the fermion velocity and the "quantum critical" aspects [10] of the interacting Dirac fermion 
problem. 

Let us begin with a statement of the problem. After taking the continuum limit to N — A 
species of two-component Dirac fermions "ifa — i ■ ■ ■ N) we have the theory defined by 
the Euclidean partition function 
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The functional integral is over fields defined in two spatial dimensions r = (x, y) and imagi- 
nary time r, cr^'^ are Pauli matrices acting on the Dirac space, and v is the Fermi velocity. 
The scalar potential which mediates the e^/(47reo|r|) Coulomb interaction between the elec- 
trons is iAT-{r,T); after a spatial Fourier transform to two-dimensional momenta q, this 
interaction requires the 2q (= 2|q|) co-efficient of the term quadratic in Ar, with the cou- 
pling g"^ = e^/tQ. The screening due to a substrate of dielectric constant e can also be 
included by modifying the coupling to [10] = 2e^/(eo(l + s))- The action S therefore 
represents the physics of an ideal graphene layer. The influence of an impurity of net charge 
Z at r = is described by 5imp- 

Many essential aspects of the theory above follow from its properties under the renormal- 
ization group (RG) transformation under which r — > r/s and r ^ r/s. A standard analysis 
shows that all three couplings in Z, namely v, Z, and g, are invariant under this transfor- 
mation at tree level. Indeed, for two of the couplings, this invariance extends to all orders 
in perturbation theory: the coupling g does not renormalize because of the non-analytic q 
co-efficient, while Z remains invariant because it is protected by gauge invariance [5]. So 
we need only examine the RG flow of a single coupling, the velocity v. Because v is a bulk 
coupling, its flow cannot be influenced in the thermodynamic limit by a single impurity, and 
so can be computed in the absence of the impurity. Such a RG flow was initially examined 
in the more general context of theories with Chern-Simons couplings in Ref. 9, but a com- 
plete presentation was given in the present context in Ref. 10: we shall use the notation 
and results of the latter paper here, with the exception that we use two-component Dirac 
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fermions with N = 4 while Ref. 10 uses four-component Dirac fermions with N = 2. 

It will be useful for our analysis to introduce two combinations of the above couplings 
which also have engineering dimension zero, and hence are pure numbers. These are 

, g'N g'Z 
^=32^ ' "=4^^ 

(we have set h — 1 elsewhere in the paper). As we will see, the coupling A is a measure of 
the strength of the electron-electron Coulomb interactions, while a measures the strength 
of the electron-impurity Coulomb interaction. 

We shall limit our explicit results here to the spatial form of the charge density 

n(r) = -5^Tr(^l(r,r)^„(r,r)), (3) 

a 

(where Tr acts on the Dirac space) induced by the impurity. However, our RG strategy can 
be extended to other observables of experimental interest, such as the local density of states. 

As noted above, we will begin in Section II by considering only the electron-impurity 
Coulomb interaction, while electron-electron Coulomb interactions will be accounted for in 
Section III. 



II. NON-INTERACTING ELECTRONS 

This section will ignore the electron-electron Coulomb interactions. Formally, we work in 
the limit A ^ 0, but a is kept fixed. The problem reduces to that of a single Dirac electron 
in the attractive impurity potential 

V(r) ^ -g. (4) 

This problem was originally studied in Ref. 4. However, they introduced an arbitrary cutoff 
at high energy to regulate the problem at short distances, and this leads to spurious results 
[5]. As we will demonstrate here, there is no dependence upon a cutoff energy scale at all 
orders in perturbation theory, provided the high energy behavior is regulated in a proper 
gauge-invariant manner. With no cutoff energy scale present, a number of results can be 
deduced by simple dimensional analysis. The Fourier transform of the charge density n(r) 
is dimensionless, and therefore we can write 

n{q) = -NF{a), (5) 

where F{a) is a universal function of the dimensionless coupling a. Note that n{q) is required 
by this dimensional argument to be g-independent, and so n(r) oc (5^(r). 

The arguments so far are perturbative, but non-perturbative effects can be deduced by 
solving the full Dirac equation in the potential in Eq. (4). This solution has appeared 
elsewhere [6-8], and so we will not reproduce it here. Such an analysis shows that the 
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FIG. 1: Feynman diagrams for the charge density without electron-electron interactions to order 
a^. The filled square is the impurity site, the wavy line is the propagator, the line is the fermion 
propagator, and the filled circle is the charge density operator. 

perturbative arguments apply for a < 1/2, but new physics appears for a > 1/2. In 
particular, Shytov et al. [6] showed that n[r) ~ — for a > 1/2 (the sign of this tail is 
opposite to that of the impurity). 

We shall limit our discussion in this section to the a < 1/2 case. One reason for doing 
so is that electron-electron Coulomb interactions act to reduce the effective value of a. This 
will become clearer in Section III, but we note here that a standard RPA screening of the 
potential V{r) in Eq. (4) can be simply accounted for by applying the mapping 

to the results of the present section. The value of A in graphene is not small [10]. 

We shall now establish the existence of the universal function F{q) in Eq. (5) to all orders 
in a. The existence of a universal F{q) is a consequence of the non-renormalization of the 
impurity charge Z [5]. We compute n{q) diagrammatically, and the needed diagrams all 
have one fermion loop and are shown in Fig. 1. 

To first order in a we have 

n{q) = -^no(g), (7) 
where Hoiq) is the bare polarization operator 

g^Nq 



16v 



(8) 



and so we have F{a) = (7r/8)a + 

The order graph in Fig. 1 vanishes by Furry's theorem, and at order we write the 
contribution to n{q) in the form 



ATf7 2^3 / d^kid^k2d^k3A(ki,k2,k3) , i , i , ^ fn\ 
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where 

A{ki, k2, ^^)=^ [{-lu + vp-a) {-tu + v{p + q) ■ a) 

X {-iu + v{p + q + ki) • a)~^ {-iu + ^(p + q + ki + ka) • a)~^] , (10) 

where it is understood here and below that — q = ki + k2 + ka. We now want to symmetrize 
this by placing the external vertex with momentum q at different points on the loop — this 
should not change the final result for n{q). In this manner we obtain 

3A{k^MM)^ J^l ^i^Tr[{-iu; + vp-a)-'{-iu; + v{p + q)-a)-' 

X {-iu + v{p + q + ki) ■ a)~^ {-iu + v{p + q + ki + k2) • 

+ Tr [{—iuj + vp ■ a)~^ {—ion + v{p + ki) ■ (x)~^ 
X {-iuj + ^(p + q + ki) • ffy^ {-iuJ + ';;(p + q + ki + k2) • 

+ Tr [{—iio + vp ■ a)~^ {—iuj + v{p + ki) • a)~^ 

X (-ia; + w(p + ki + k2)-a)"^(-ia; + u(p + q + ki + k2)-a)"^]|. (11) 

Now this expression has the important property that it vanishes at q = 0, where we have 

3A(ki,k2,k3) = /0/ [{-iu-rvp-a)-'{-iu-rv{p-r\i,)-a)-^ 

X (-ia; + 'i;(p + ki+k2) -a)"^] . (12) 

This property allows us to establish that the integral in Eq. (9) is convergent and cut- 
off independent. Let the loop momenta /ci, A;2, and ^3 all become much larger than the 
external momentum q. The resulting integrand will scale as the power of momenta associated 
with a logarithmic dependence on the upper cutoff. However, in this limit of small q we have 
just established that the integrand is zero. It is clear that this argument can be extended to 
all orders in a. We have thus established the existence of the cut-off independent function 
F{a). We computed the integral in Eq. (9) numerically, and so obtained 

F{a) = + (0.19 ± O.Ol)^^ + O(a^). (13) 



III. INTERACTING ELECTRONS 



We will now consider the full problem defined in Eq. (1), and account for both the 
electron-electron and electron-impurity Coulomb interactions. 

The problem can be solved in two limits: in the weak coupling limit A — > and the large 
N limit, N ^ 00 with fixed Z = 0(1). In both cases a/{l + X) <^ 1, so one can hmit oneself 
to linear response in which the induced charge is [generalizing Eq. (7)] 



n{q) = -ZD{q)U{q), 



(14) 
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where D{q) is the full propagator of the Coulomb potential At, and Il{q) is the polarization 
tensor. The connection between D{q) and n(g) is 

D-\q) = D^\q) + U{q), (15) 

where Do{q) is the bare propagator, 

^o(g) = ^. (16) 

To leading order (either in couphng or l/N), the polarization operator was given in Eq. (8), 
and we showed in Section II that this gives rise to a independent n{q), or a screening 
charge localized at r = 0. 

However, if we compute corrections, wc find logarithmically divergent diagrams, where the 
logarithms are cut off from above by the inverse lattice size and from below by q. The leading 
logarithms are summed by a standard RG procedure. Since the theory is renormalizable, 
we can eliminate the dependence on the cutoff by expressing the each diagram in terms of 
the renormalized parameters, instead of the bare parameters of the Lagrangian. Choosing 
the renormalization point to be qo, and denote vo as the fermion velocity at the scale v, the 
polarization tensor can be schematically written as 

U{q)^U{q;qo,Vo). (17) 

In n there are logarithms of the ratio q/qo- We notice that Il{q; qo, fo) is is invariant under 
a change of the renormalization q^,, given that Vq is changed correspondingly (the particle 
density has no anomalous dimension). To eliminate the powers of log{q/qo) we can choose 
qo — q, hence 

U{q)^U{q;q,v{q)), (18) 

where in the perturbative expansion of the right hand side there is no large logarithms. Thus 
to leading order it is given by a single diagram, which was computed previously [Eq. (8)], 

All the leadings logarithms are contained in the function v{q), which satisfies the equation 

q§^v{q)^P{v), (20) 

with the boundary condition v{qo) — Vq. The screening charge is then 

niq) = -Z-^ , A(g) = ^ . (21) 

The problem is now reduced to the problem of finding v{q) [or, equivalently, A(g)]. This 
problem has a long history [11]; most recently it has been revisited in Ref. 10 (see also 
below) . 
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To find the spatial charge distribution n(r) one needs to take Fourier transform of Eq. (21). 
First one notice that if the velocity does not run then n(r) is proportional to S{r). Only 

when V runs with the momentum scale does n{r) differ from away from the origin. When 
the running is slow (as at weak coupling or at large A^), the amount of screening charge 
enclosed inside a circle of radius r (assumed to be much larger than the lattice spacing), to 
leading order, is 

r 

A(g) 



J dr'n{r') ^ n{q)\q=i/r = 



(22) 

q=l/r 



1 + m 

The total screening charge is small if A at the scale 1/r is small, and close to —1 if A is large. 
Differentiating both sides of Eq. (22) with respect to r, one finds 

nir) = ^(^) (M^ (23) 

^ ' 27rr2 [1 + A(g)]2 v{q) ' ^ ' 

Note that the beta function for v is negative, therefore we arrive to a counterintuitive result 
the screening charge is positive. To see what is happening, let us take the limit r — > oo in 
Eq.(22). This hmit corresponds to the infrared hmit — > 0. We know that asymptotically 
v{q) grows to oo in this limit (although only logarithmically), hence 



oo 

/ 



dr'n(r') = 0. (24) 



i.e., the total screening charge is zero when integrated over the whole space (although the 
integral goes to zero very slowly). The presence of an external ion, therefore, only leads to 
charge redistribution: a fraction of the unit charge is pushed from short distance (of order 
of lattice spacing) to longer distances, but none of the charge goes to infinity. Therefore, 
there is a finite negative screening charge localized near r = 0. Its value can be found by 
taking r to be of order of inverse lattice spacing in Eq. (22). The final result for the 
screening charge density can be written as 

^iTW^r"^ 2^ [1 + \{q)f ^JioT ■ ^^^^ 

In the rest of the note we will concentrate our attention on the long-distance tail of n{r), 
ignoring the delta function at the origin. 

At weak coupling (A -C 1), the beta function for v{q) is 

= (26) 
The solution to the RG equation, with the boundary condition v — Vq at q — qo, is 

v(q) = + ^ In ^ , (27) 
Idtt q 
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and the screening charge density is 

Notice that the resuh is proportional to the square of the small coupling constant A = 
g'^N/32v, although we have performed the calculation to leading order in the coupling. The 
reason is that for the charge density n{r) to be nonzero, it is necessary that the coupling 
constant runs. The density n(r) therefore contains the beta function P{v), as seen in Eq. (23), 
and hence is second order in the coupling constant. 

In the 1/N expansion the beta function for v{q) was computed in Ref. 10: 



8. /ln(A + ^A^)^^ . ^^^^ 



P{v) = <^ 



( arccos A tt , 

+ 1-:tt , A< 1. 



TT^iV VA^T^ 2A 



(29) 



The two expressions smoothly match each other at A = 1. 

In is instructive to analyze two regimes where the RG equation can be solved analytically. 
The first regime is A ^ 1 where the result is the same as in Eq. (28). The second regime 
is the strong-coupling regime A ^ 1. This regime corresponds to a quantum critical point 
characterized by a dynamic critical exponent whose value at large N is [10] 

. = l-^ + 0(iV-). (30) 

In this regime (5 — {z — l)v. The solution to the RG equation, with the initial condition 
V — Vq aX, q — Qo, is 

„(,) = „„(^)""', l-z^^. (31) 

In this regime 

i.e., the charge density follows a power law behavior n(r) ~ r'^~^. The power is slightly 
different from —2. 

In real graphene A is of order 1, so one has to solve numerically the RG equation. We 
chose the scale Qq to be comparable to the inverse lattice spacing, r^^, and Vo to be lO^m/s, 
a typical value found in experiments. We then run v according to the leading (in 1/A^) RG 
equation in two cases, in vacuum and when graphene is on a Si02 substrate with dielectric 
constant e = 4.5. We then plot 27rr^n(r) as a function of the distance r on Figs. (2) and (3). 

As seen from the figures, the charge density n(r) roughly follows the r^^ law: when r 
changes by two orders of magnitude, the product r'^n{r) changes by a factor of less than 1.5 
in both cases. 
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logio(r/ro) 

FIG. 2: The dependence of 27rZ~^r^n(r) on the distance r for suspended graphene. Note that 
coordinate r is on a logarithmic scale. 




Iogio(r/ro) 

FIG. 3: The dependence of 27rZ~^r'^n{r) on the distance r for graphene on a substrate with e = 4.5. 
Note that coordinate r is on a logarithmic scale. 

IV. CONCLUSIONS 

In this paper we have considered the problem of screening of a Coulomb impurity 
in graphene. We show that there is a qualitative difference between screening by non- 
interacting and interacting electrons. In the case of non-interacting electrons the induced 
charge density is localized at the position of the impurity when the impurity charge is small. 
The interaction between electrons lead to a long-distance tail in the induced charge distri- 
bution, with a counterintuitive sign which is the same as that of the impurity. 

One problem that is not addressed in this paper is the screening of an impurity with 
large a ~ 1 by an interacting electron gas. We hope to address this problem in a future 
publication. 
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